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Abstract. This is a companion report for the paper "Transverse electric scattering on inhomo- 
geneous objects: Spectrum of integral operator and preconditioning" by the present authors [6]. In 
this report we formulate the two-dimensional transverse electric scattering problem as a standard 
singular integral equation, derive the symbol of the integral operator for Holder-continuous contrasts, 
and calculate the elements of the system matrix obtained after discretization via the mid-point rule. 
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1. Derivation of the TE singular integral equation. Although some of 
the material presented here may be new, in particular, the explicit expression for 
the operator symbol, the methods and techniques are standard and can be found in 

mum- 

Starting from the frequency-domain Maxwell's equations [B] with the time con- 
vention e~ tult we arrive at the following set of equations for the scattered fields with 
the induced currents as sources: 
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where 



J^ d (x, u) = -iw [e(x, u) - £ b ] B fc (x, u), k = 1, 2, 3; 
i4T™ d (x,w) = -iu [/i(x,w) - fi-b] H m (x,uj), m = 1,2,3. 



(1.3) 



Equation (JTTTJ) describes the TE case while (11.2[) the TM case and is the dual of (|l.lj) . 
The two-dimensional Fourier transform of Ijl.ljl and ()1.2|) with respect to coordinates 
x\ and x 2 takes us from the (x, uS) domain to the (k, uj) domain, whereby d n — > 
—ik n , n — 1, 2. In this way we arrive at the following linear algebraic problem of the 
form AF = — S for the TE case: 
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Solving (jl.4[) via matrix inversion, we obtain the scattered fields in terms of the 
induced currents: 
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where we have introduced the vector potentials 
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m = 1,2,3. 
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Transforming (11.51) back to the (x, w) domain and recognizing the partial derivatives 
as k n — ¥ id n , n — 1,2, we get 
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The (x, w)-domain vector potentials are the spatial convolutions of the induced cur- 
rents with the scalar Green's function, namely 



A k (x, uj) 



x'ei 



x'ei 



3 (x-x', W )4" d (x, W ), fc = l,2,3; 
9 (x-x', U )C d (x, W ), in — 1,2,3. 



(1.8) 



The Green's function g(x, u>) is the two-dimensional inverse Fourier transform of the 
expression l/(fc 2 +fc§ — w 2 £b/^b) appearing in (|1.6I) . It is easy to verify that this Green's 
function satisfies the non-homogeneous Helmholtz equation with a line current (two- 
dimensional Dirac's delta function) as a source term, located at x' position. Since 
the scattered fields are supposed to satisfy the so-called radiation boundary condition 
(i.e. outgoing waves decaying at infinity), out of the two possible solutions of the said 
Helmholtz equation one choses 



s(x-x>) = |ff«(fc b |x-x'|). 



(1.9) 



The other possible solution has the form of the Hankel function of the second kind, 
and is chosen when a different time convention is used, i.e. for the time-dependence 
of the form e ZUJt . 

Substituting the induced currents from (|1.3[) and expressing the scattered fields 
as E sc = E — E m and H sc = H — H ln , we arrive at the following integro-differential 
equations with the total fields as the fundamental unknown: 
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where 

Xe (x,w =-^-2-^-1, (1.11) 
Xm (x >W ) = ^-l J (1.12) 

Mb 

are the normalized electric and magnetic contrast functions, respectively, and the star 
(*) denotes the 2D convolution, for example 

g*UcEi)= I <?(x-x',u;)xe(x')£i(x')crx'. (1.13) 

Jx'GK 2 

The system (jl.101) contains nine scalar integro-differential operators. If the partial 
derivatives are carried out, then we arrive at nine "pure" integral operators whose 
kernels have different degrees of singularity. The weakly singular kernels result in 
compact operators, whereas the strongly singular kernels require some extra caution. 
The presence of the second-order partial derivatives in the upper left (2 x 2)-corner 
of the derivative matrix in (| 1 . 10|) indicates that we have four scalar strongly singular 
kernels in the TE case. The corresponding operators are called singular integral 
operators. 

To arrive at the standard form of the singular integral operator we first separate 
the domain of integration D in two sub-domains as D = [D \ D(e)] U D(e) where 
D(e) is a circular area around x with the radius e. Let u r be either of the electric 
field components E\ or E2, then the product between the operator matrix and the 
convolution vector in (|1.10j) can be written as 



h = lim / d k d r g{x - x')xc(x')ur(x')dx' 

+ limd k [ <9 r . 9 (x-x')xc(x>r(x')dx', k,r = l,2. (1.14) 

Here and in what follows we save some space by not showing the parametric depen- 
dence on uj. The first term in (|1.14j) is recognized as the principal value, while the 
second integral will be denoted by li- This second term incorporates the Green's 
function of (|1.9j) . Utilizing asymptotic expansions for small arguments for the zero 
order Bessel and Neumann functions PQ, namely 



J (z) = l-j + O(z% (1.15) 
N (z) = -jJ (z) + - J (z)ln 7 - + 0(z 2 ), (1.16) 

7T 7T 2 

we replace the Hankel function in Green's function by its asymptotic expansion based 
on (j!.15[) and (|1.16p . arriving at 

is 7 i 1 fcyjx — x'l ,, 0> 
5 (x-x0^ + i --ln^^ + O(|x-xf) 

= jo(x-x')+Si(x-x'), (1.17) 

where we have defined 



.go(x - x') = g(x - x') - 



1 , fcb|x-x'| 
— —In 
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and 

1 , &hlx — x'l 

9l (x-x') = -— ln-^i— L (1.19) 

The function go(x — x') is not singular and its contribution (even after differentiation) 
in I2 will vanish in the limit e — > 0. The nonzero contribution to J2 comes from 
<7i(x — x'). Taking the first partial derivative d r of gi(x — x'), we get 

^(x-x')^^, (1.20) 

where O r = (x r — xj.)/|x — x'|. A singularity of order one has appeared, which is 
still a weak singularity (the order of the singularity is smaller than the dimension of 
the manifold, which is two-dimensional in the present problem). Applying the second 
partial derivative dk we obtain a strong second-order singularity: 

d k d r9l (x - x') = ^0 k ^ = ±.^L^ (29,6,. - S kr ) . (1.21) 

Thus, omitting the go part (since it disappears in the limit) we re-write the second 
term of (fTTI)) as 



\imd k / 9 r g(x - x')Xc(x')ur(x')dx' 

= lim f ^ — %Xc(x'K(x')dx'. (1.22) 

^ix'PDW 27r l x - x l 



c'ei>(£) 

By adding and subtracting x c (x)w r (x) inside the integral of (|1.22[) . we get 



Mmdk 9 r ^(x - x')xe(x')Mr(x')dx' 

Jx'eD(e) 
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1 f — 9 r , 
lim/ d k - —t [xe(x')ii r (x') - Xe(x)u r (x)] dx' 

' " Ve£>(e) l x - x I 



27T £->0 



+ ^lim X c( x K( x ) / d k — %dx'. (1.23) 

27T^0 J*><=D(e) l x - x l 

Assuming the Holder continuity of the function x c (x)it r (x), i.e., assuming that there 
exist a, C > 0, such that for all x, x' el 2 , 

|Xe(xR(x) - X c( x 'K( x ')l < C|x - x'| Q , (1-24) 

we effectively lower the order of singularity in the first integral in the RHS of (|1.23j) . 
Hence, in the limit e — > this term is zero. Interchanging dk with —d' k and applying 
the 2D divergence theorem in the remaining term of (ll.23[) , we obtain 

\imdk / <9rff(x - x')xe(x')ttr(x')dx' 

J x >eD{c) 

1 / -0 r 
= ^Xe(x)u r (x) Km i> -v k - ^dx', (1.25) 

2?T e->0 J x / e9£ ,( e ) |X-X'| 
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with v k = — Gfc being the projection of the normal unit vector of the polar coordinate 
system on the axis. Finally, we arrive at the result 

limdk / d r g(x — x')xe(x')it r (x')dx / 
J X 'eD{c) 

= ^ Xc (x)^(x)lim£^ r ^dx' 

= -7r~xo(x)u r (x) iim * e fc e r dx' 

= -ixc(x)ur(x)5 rfc , (1-26) 



where we notice that the last contour integral is over the unit circle. Hence, (|1.14[) 
can now be written as 

h=p.v. I d k d r g{x- x')Xo(x')ur(x')dx' 
Jx'eD 



~ 2^ c W u, '( x ^ rfe ' k,r = 1,2. 

Now, we can represent ()1.10p in the standard form 
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(1.27) 







~S 





0" 








Ef 







s 







E2 


+ p. V. 












I 




#3. 





X e Ei 
X C E 2 
X m H 3 

X c Ei 
X C E2 
X m H 3 



(1.28) 



with 5 denoting the operator of (pointwise) multiplication with the function s(x) = 
l + l/2x c (x), and / - the identity operator. The kernels in the principal-value operator 
in (|1.28[) are easily recognized using (| 1 . 2 1 1) and (| 1 . 10[) . and are given by 

G„ m (x) = - r [28„9 m - S nm ] , n,m = 1, 2. (1-29) 
Z7r|xp 

The kernels in the compact operator in (|1.28[) are easily obtained using (|3.3[) . (|3.4p 
and (ll.lOp . and are given by 
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2. Derivation of the Symbol. Using operator notation, (|1.28[) can be written 



(AE)(x) = ( E + -M ) E(x) + A s ME(x') + KE(x'), 



(2.1) 



where A s is a (matrix) singular integral operator, K is a (matrix) compact integral 
operator, E is a (matrix) identity operator, and M(x) = Xc(x)E2 is a (matrix) mul- 
tiplication operator, where the lower-right element of E2 is set to zero. The singular 
term of (12.11) is of the form 



(i s ME)(x) 
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According to (|1.21[) . the characteristic (matrix) function F(0) has the form 

1 



F(0) = - 



2tt 
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(2.3) 



where I2 is the (3 x 3) identity matrix with the lower-right element set to zero, and 
the tensor Q is given by 
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(2.4) 



For Xc(x) Holder-continuous on R 2 , the symbol of the compound singular integral 
operator of (|2.1j) can be computed as (see [3]): 



Smb(A) = I + Jxe(x)I 2 + Smb(A s ) X c(x), 



(2.5) 



where I is the ordinary (3 x 3) identity matrix. Since the characteristics F depends 
only on x — x', the symbol of the singular integral operator A s is the Fourier transform 
of its kernel F(0)/|x— x'| 2 with respect to the variable y = x — x', i.e. it is a k-domain 
matrix- valued function F s (k). Computing this Fourier transform is a daunting task, 
and we shall use a shortcut proposed in [4]. 

Let A s be a single component of our matrix- valued operator A s and let $ s (k) de- 
note its symbol, which is one of the components of the matrix- valued symbol function 
F s (k) we are trying to compute. The symbol Smb(A s ) = $ s (k) of a scalar simple 
singular integral operator can be expanded in a series of 2D spherical functions of 
order p, that is, in Fourier series of sines and cosines [3], namely 



* s (0) 



p=0 



(2.6) 



where = k/|k| is the unit vector in the k-domain, — sm(p4>), and Y^lj 

cos(p4>) is the basis of the expansion, and are the expansion coefficients, and 
4> is the directional angle of the unit vector 0. In (|2.6p . 72,p = 7ri p T(p/2)/T((2+p)/2) 
[3]. Since any component /(0) of the characteristic matrix- valued function F(0) 



-(2) 
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given by 
well 



depend only on 0, we can expand each of them in a Fourier series as 



oo 

/(©) = [< } sin(W>) + cos( P 0) 



(2.7) 



P =i 



where is the directional angle of the unit vector = (x — x')/|x — x'|. It was shown 
in [J] that the expansion coefficients in (|2.7I) are the same with those of (|2.6p . Using 
(|2.3p and (|2.4I) . we see that the components of the characteristic F are 



F(0) 
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'2tt 



x - x') - h] = - 



2lT 
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(2- 



To calculate the expansion coefficients, we must equate each element [F(0)]fc r , k,r = 
1,2, from (|2.8p with the series from (|2.7[) . For example, for the element [F(0)]n we 
have 



1 

27 



[2 cos 2 



1 

2^ 



OO 

cos(20) = ^ sin(p0) + a(, 2) cos(p0) . (2.9) 



P =i 



From this expansion it is obvious that the only nonzero expansion coefficient is 



,(2) 



T/(27r), while the rest are all zero (i.e. 



Vp and 4 2) = Vp ^ 2). The 



same procedure is followed for the rest of the components in (|2.8I) . Then, we substitute 
the known expansion coefficients in (|2.6[) to get the elements [F s (0)]/ cr , k, r = 1,2. 
Continuing our example, the element [F s (0)]n is obtained through the substitu- 



tions 72,2 = 7Ti 2 r(l)/r(2) = -7T, Y 



(2) 



cos (20) and a. 



(2) 



T/(27r), and therefore 



[F s (0)] n = l/2cos(20) 
ponents, we finally get 



= F S (0) = 



cos 2 (0) — 1/2. Following the same procedure for all com- 



cos 0—1/2 sin cos 
sin cos sin 2 0—1/2 




= Q(k) - -i 2 . 



(2.10) 



Substituting the result of (|2.10[) back in (|2.5[) . we get the symbol of the complete 
operator as the following (3 x 3) matrix-valued function: 



Smb(A)(x,k) =I + Xe (x)Q(k). 



(2.11) 



3. Derivation of the algebraic system matrix. Equation (|1.28j) defines an 
algebraic system Au = b for the numerical evaluation of the fields. To derive the 
matrix elements, its more convenient to use the equivalent integral form of (jl.28l) . 

The application of the matrix operator of (jl.101) -the matrix that contains the 
partial derivatives- on Green's functions, give us the Green's tensor. In order to get 
the equivalent integral form of (|1.28[) . we need to split the Green's tensor in two 
parts. The first part, denoted by A(x — x'), corresponds to the second order or mixed 
derivatives only, i.e. 
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The second part, denoted by B(x — x'), corresponds to the first order derivatives only, 
i.e. 



x- x' = 



d 2 ' 
-di 
-di di 



.g(x-x'). 



(3.2) 



To obtain the explicit relations for A(x — x') and B(x — x'), we first calculate the first 
order derivative which is 



d r g(x - x') = -ifcb^-^ffW^lx - x'|), r = 1, 2. 
4 x — x' 



(3.3) 



Then, the mixed derivative is 
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where 8kr is the Kronecker's delta. So, the A tensor is given by 



ffp^fcblx-x'l) 
fc,r = l,2; (3.4) 



(x - x') =lfc b ^-H 1 (1) (fc b |x - x'|)[2Q(x - x') - l a ] 
4 |x — x'| 



-^if^(fc b |x-x'|)[Q(x-x')-I 2 ]. 



(3.5) 



The tensor Q was introduced in (|2.4[) while la is, as explained in Section^ the (3 x 3) 
identity matrix with the lower-right element set to zero. The B tensor is obtained 
with the use of (13.31) 



B(x-x') = --k h H\ 1) (fc b |x-x'|)0(x-x')x, 
where we have now introduced the 0x tensor given by 

X2 ~ x 



(3.6) 



0(x-x')x = 
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X — x' | 

X\ —x^ 
|x— X 7 I 







(3.7) 



|x-x'| |x-x'| 

After this splitting, we can easily express (|1.28[) in the following integral form 
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where the vectors E(x) = [-Bi(x), E 2 (x), 0] T and H(x) = [0, 0, i? 3 (x)] T . Same applies 
for E in (x) and H in (x). 

Equation (|3 . 8|) can be rewritten as 



E m (x) 



l + ^Xc(x) 



E(x) 

-p.v.f A(x-x')Xe(x')E(x')dx' 
- f A(x-x')Xc(x')E(x')dx' 

Jx'£D\D s 

-iwfi h [ B(x-x')Xm(x')H(x')dx' 

-iuj» b [ B(x-x')Xm(x')H(x')dx', (3.10) 

Jx'£D\D B 

where we have separated the domain D into a domain D \ D s "free" of singularity, 
and a domain D s which encloses the position vector x and hence the singularity. In 
this way we have three integrals in the usual sense and one in the sense of principal 
value. Applying a simple collocation technique with the mid-point rule, the usual 
sense integrals over the domain D \ D s would be given by 



L 



A(x-x')Xc(x')E(x')dx' 

'ED\D S 
N-l 

A(x„ - x m )xo(x m )E(x m )S* m , n = l,2,...,N; (3.11) 

m— 1 



iu)fi h [ B(x - x')Xm(x')H(x')dx' 



i'GD\D B 
N-l 

sa -iujn h 2J B(x„ - x m )xm(x m )H(x m )5 m , n = 1, 2, . . . , N; (3.12) 

m=l 

where S m — h 2 is the surface of each elementary cell D m of side h in the computational 
domain. 

We now proceed to calculate the principal value integral in p.lO|) . We have 



p.v. [ A(x-x')Xc(x')E(x')dx' 
= lim / A(x-x')Xc(x')E(x')dx' 



WeD s \D{e) 

PS Xc(x„)E(x„) lim / A(x - x')dx', (3.13) 

Jx'eD s \D{e) 

where the domain of integration in the last integral is over the singular cell but with 
an exception of a small circular neighborhood having radius e. We further define the 
tensor L as 

L = lim / A(x - x')dx'. (3.14) 



e->0 



■eD s \D(e) 
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Now we recall the form that A(x — x') has from (13.51) . It is convenient to rewrite 
in terms of Q(x — x') and I 2 , namely 



lx-x j=i 



2k h - 



-^(fcblx-x'D-fe^W^lx-x' 



+ 



klH^(k h \x - x'|) - k hr ^-H[ L) (k h \x - x'|) 



C 1 )/ 



>(x - x') 

h- (3.15) 



To compute (|3 . 14|) . we transform the singular cell D s from rectangular shape to a 
circular disk having center at x and radius a n . The transformed circular cell with the 
original one have the same surface S = ira^ = h 2 . Then, (|3.15|) is transformed into 
polar coordinates as 



1 



2k h -H i 1 1) (k h p)-kiH^>(k h p) 
P 



k 2 ,H ( a 1 \k h p)-k h ^H y 1 1 >(k hP ) 



(3.16) 



Now, we calculate each element of the tensor L from 



2k h - p H[ 1 \k hP )~klH^>(k hP ) 



<-2rr(l) 



p=e 



klH^{k h p)-k h -^\k h p) 



P&P / Qpq((t))d(j) 
J<f>=0 

pdp I 6 pg d<t>\, p,q = 1,2. (3.17) 

0=0 J 



The integration in (13.171) is carried out only three times since the tensor Q(x — x') is 
symmetric. The angular integral with integrand function Q pq ((t>) gives ir6 pq , while the 
other with the integrand function 8 pq gives 27r<5 pg . The radial integrals are evaluated 
by changing variables z = k^p and using the well known integral formulas for Bessel 
functions [5] 



f 



C" +1 Z„(C)dC = ^ +1 Zn+lM, 



r n+l z„(c)dc--z-" +l z„_ 1 (z). 



(3.18) 
(3.19) 



The result is 



£->0 



-k h H[ 1 \k h a n )~—k h lim eH[ X \k h e) \, p,q=l,2. (3.20) 



The involved limit is easily calculated giving —i2/(nkb). Therefore, substituting the 
results for the tensor L back to (|3.13[) 



p.v. f A(x-x')Xe(x')E(x , )dx' 

Jx'GD, 



Xe(x„) 



1 iira r 



-khH^ '(Vn) 



E(x„), x„ e D s 



(3.21) 
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What remains is the calculation of the integral involving B(x— x') over the domain 
D s in (03) ■ We have 

-iujfih / B(x - x')x m (x')H(x')dx' 

« -^ b x m (x n )H(x„) f B(x-x')dx'. (3.22) 

Utilizing (|3.6p and (|3.7[) . and expressing B(x — x') in polar coordinates, in the same 
manner as we did for the L tensor before, we arrive at an angular integration from 
to 2ir for each component Q pq of 0(x — x') x tensor. This integration is zero for each 
component O pg , and hence this term does not contribute. 

Putting together the results of (|3 . 1 ljl . (|3.12[) and (|3.21[) . we arrive at the linear 
equations 

jl - Xc(x„) 
N-l 

- h 2 ^ A ( x " _ x m )xc(x m )E(x r „) 

in— 1 

N-l 

- iujfi h h 2 B ( x - - x m ) Xm (x m )H(x m ) = E in (x n ), n = 1, 2, . . . , N. (3.23) 

m— 1 

Following the same analysis for (I3.9[) . we arrive at the linear equations 

jl - Xm(x„) 
N-l 

- h 2 ^2 M*n - X m )Xm(x m )H(x m ) 

m=l 

N-l 

+ iue h h 2 B(x„-x m ) X o(x m )E(x m ) = H in (x„), n = 1, 2, . . . , N. (3.24) 

m— 1 

Equations f|3.23D and (|3.24[) compose the linear system for the evaluation of the un- 
known total fields, with the following structure 



An 


A 12 


Ais 




Ul 




~h~ 




An 


A 2 2 


A 23 




U2 




b 2 


(3.25) 


A 3 i 


A 32 


A 33 




"3 




b 3 





In (|3.25[) . [m, U2, u 3 ] T = [Ei(x n ), 2£ 2 (x„), iJ 3 (x„)] T , n — 1, 2, . . . , N, are the unknown 
total field components on the grid, while [bi,b 2 ,b 3 ] T = {E^ (x n ) , E™ (x n ) : H^ 1 (x n )] T 
contains the grid values of the incident field. The elements of the system matrix in 
(|3~2"5j) are now easily recognized with the use of (pHHl) . (PT2"4]l . (|53j) . (|31)1) and 



iirk^h 
40F 



E(x„) 



ink^h 
2^ 



H(x„) 
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The results are 

[A U ] nm = -fcb^ 2 Xe(x TO ) 



Z/C b r nm 4 



(^b^nm) 



5 u y 1=1,2, m^n; (3.26) 



[Au] nn —1 



nrk h h (1) _ 
4V7T 



Xe (x n ), *=1,2; (3.27) 



[Acqlnm — — k^ti 2 Xc(~X-rn) 77 -^1 '(^b J 'tmi) ~ "7^0 ^ (^b r nm) 

X ^,nm^,Tim(l -<5nm), £,q=l,2, I ^ q; 



(3.28) 



[-4., 



iw£b^ 2 Xo(x m )^i/ 1 (1) (fc b r rim )6'i,„ m (l - 5 nm ); 



[A 3 l]nm =i^£bh 2 X C {^rn)^H[ 1 \k h r nm )9 2 , nm (l - Sum)] 



[A 2 z\ nm = - «w/i b /i 2 Xm(x„ l )^ff 1 (1) (fc b r„ m )6 l i,„ m (l - S nm ); 



(3.29) 
(3.30) 
(3.31) 
(3.32) 



(3.33) 



[A 33 }nn = 1 + 



1 - ffj; >{k h h/y/ir) 

Zy'TT 



Xm(x„). 



(3.34) 



In (|3.26l) -p.34p we have denned r nm = |x n -x m |, 9 Lnm = (a^,„ -a^, OT )/r„ m , 1= 1,2, 
n, m = 1, . . . , N, and xi >n denotes the Cartesian component of the 2D position vector 
x n , pointing at the nth node of the grid. 
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